'data.frame': 26 obs. of 3 variables:
$ seedlings: int 7 10 0 5 54 0 30 14 1 2 ...
$ dungs : int 5 0 330 960 115 2330 475 10 2000 130 ...
$ property : Factor w/ 2 levels "Campesino","Empresa Forestal": 2 2 2 2 2 2 2 2 2 1 ...
Hora
9:30 a 11:00 teoría
11:00 a 11:30 ejercicios
11:30 a 11:45 break
11:45 a 12:45 teoría
12:45 a 13:30 ejercicios y dudas
Temario
GLM
Ejercicio GLM
Análisis multivariantes
Ejercicio PCA + GLM
Definición
Estructura
2.1. Función de vínculo
2.2. Familia de distribución de errores
Supuestos modelos lineales:
Los errores se distribuyen normalmente
La varianza (residual) es constante
La variable respuesta se relaciona linealmente con la(s) variable(s) independiente(s), si alguna de estas variables es continua
Uno o varios supuestos no se cumplen.
¿Qué hacer?
GLMs: extensión modelos lineales especificando estructuras no normales de errores y varianzas no constantes.
Modelos lineales caso particular de GLMs con distribución de errores normal y varianza constante.
Tipos de variables que por definición violan supuestos modelos lineales:
conteo de casos (p.ej. abundancia de especies)
conteo de casos como proporciones (p.ej. porcentaje plántulas muertas en un experimento)
respuesta binaria (p.ej. usado o no usado, presente o ausente)
3 componentes principales:
1. El predictor lineal
2. La función de vínculo
3. La estructura de los errores
Estructura modelos lineales:
yi = β0 + β1x1 + … + β1xn + Ɛi
Ɛi ~N(0,σ2)
Estructura GLMs:
f(yi) = β0 + β1x1 + … + β1xn + Ɛi
Var(y) ~ f(µ)
Transformación de la variable y buscando linealizar relación y con x y/o transformando en continuas variables que no lo son.
Por tanto, cambia la naturaleza de la variable respuesta.
Determinadas funciones de vínculo deben implementarse con determinados tipos de variables respuesta
Poisson: varianza aumenta linealmente con la media. Para respuestas tipo conteo.
Binomial negativa: extensión de Poisson. Para variables conteo con sobredispersión.
Binomial: para proporciones y datos presencia/ausencia.
Gamma: para datos continuos con coeficiente de variación constante.
Función glm()
Argumento family:
glm(y ~ x, data = data, family = “family”)
Datos de un estudio sobre el efecto de la ganadería y dos tipos de ganadería sobre la regeneración de la araucaria. En este trabajo se muestreó el número de plántulas de araucaria en 36 parcelas repartidas entre dos tipos de ganadería: pequeños propietarios, que hacían un uso más intensivo de la tierra, y grandes empresas forestales, que mantenían dentro de sus explotaciones algunos remanentes de bosque nativo. La actividad del ganado se midió de forma indirecta contando el número de bostas en cada una de las parcelas.
Variables:
seedlings: Número de plántulas de araucaria
dungs: Número de bostas (proxy actividad del ganado)
property: Dos tipos de ganadería: pequeños propietarios y grandes empresas
Hipótesis:
H0: β = 0 (pendiente entre plántulas y número de bostas es 0)
H0: µ1 = µ2 (número plántulas en pequeños propietarios y empresas sin diferencias)
H0: β1 = β2 (relación entre plántulas y bostas es igual en pequeños propietarios y empresas)
glm_arau1 <- glm(seedlings ~ dungs*property, ara, family = "poisson")
anova(glm_arau1, test = "Chi")Analysis of Deviance Table
Model: poisson, link: log
Response: seedlings
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 25 347.56
dungs 1 30.967 24 316.59 2.624e-08 ***
property 1 80.598 23 236.00 < 2.2e-16 ***
dungs:property 1 37.974 22 198.02 7.170e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Devianza explicada por el modelo
Medida de cantidad de variabilidad que explica el modelo (similar R2 en modelos lineales)
D2 = 1 - Dmodelo/Dnulo
Call:
glm(formula = seedlings ~ dungs * property, family = "poisson",
data = ara)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.497423 0.135408 18.444 < 2e-16 ***
dungs -0.010211 0.001789 -5.706 1.15e-08 ***
propertyEmpresa Forestal 0.576261 0.170467 3.380 0.000724 ***
dungs:propertyEmpresa Forestal 0.009041 0.001803 5.015 5.31e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 347.56 on 25 degrees of freedom
Residual deviance: 198.02 on 22 degrees of freedom
AIC: 282.3
Number of Fisher Scoring iterations: 5
Fórmula: log(y) = 2.497 + 0.576 Empresa + (-0.010 + 0.009 Empresa) Bostas
Call:
glm.nb(formula = seedlings ~ dungs * property, data = ara, init.theta = 1.13715532,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.152593 0.362295 5.942 2.82e-09 ***
dungs -0.005663 0.002468 -2.294 0.0218 *
propertyEmpresa Forestal 1.113371 0.557993 1.995 0.0460 *
dungs:propertyEmpresa Forestal 0.003945 0.002528 1.561 0.1186
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(1.1372) family taken to be 1)
Null deviance: 49.506 on 25 degrees of freedom
Residual deviance: 28.040 on 22 degrees of freedom
AIC: 155.22
Number of Fisher Scoring iterations: 1
Theta: 1.137
Std. Err.: 0.384
2 x log-likelihood: -145.218
2 variantes de estos modelos:
Variable con valores binarios (1/0): modelos de regresión logística
Variables con valores porcentuales (nº éxitos / nº de ensayos)
Ejemplos:
Función de vínculo:
logit: log(p/(1-p))
logit(y) = ln(p/(1-p)) = β0 + β1x1 + … + β1xn
Curva sigmoidea entre 0-1
Datos de un estudio que investiga la estrategia reproductora de una especie invasora de alga en relación a dos especies nativas. Una de las preguntas es si la probabilidad de encontrar estructuras reproductoras para individuos de la misma talla es diferente entre la especie invasora y las dos especies nativas.
¿Probabilidad de encontrar estructuras reproductoras es diferente entre la especie invasora y las nativas teniendo en cuenta la talla del individuo?
Hipótesis: a misma talla del individuo, la especie invasora tendrá mayor probabilidad de reproducirse que las nativas
Variable respuesta (Estado):
- 0: sin estructuras reproductoras
- 1: con estructuras reproductoras
Variables predictoras:
- talla del individuo (cm)
- especie (Sp)
Analysis of Deviance Table (Type III tests)
Response: Estado
LR Chisq Df Pr(>Chisq)
Long 21.1448 1 4.259e-06 ***
Sp 3.6498 2 0.1612
Long:Sp 0.4590 2 0.7949
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Interacción no significativa
Analysis of Deviance Table (Type III tests)
Response: Estado
LR Chisq Df Pr(>Chisq)
Long 100.112 1 < 2.2e-16 ***
Sp 31.771 2 1.262e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] 0.2829254
Call:
glm(formula = Estado ~ Long + Sp, family = "binomial", data = algas)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.80703 0.42701 -6.574 4.91e-11 ***
Long 0.26505 0.03174 8.351 < 2e-16 ***
SpTom -0.57476 0.34528 -1.665 0.096 .
SpVer -1.84877 0.36933 -5.006 5.57e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 476.38 on 358 degrees of freedom
Residual deviance: 341.60 on 355 degrees of freedom
AIC: 349.6
Number of Fisher Scoring iterations: 5
x_seq <- seq(min(algas$Long), max(algas$Long), length.out = 100)
newdata <- expand.grid(Long = x_seq, Sp = levels(algas$Sp))
# Obtener predicciones en escala de respuesta (probabilidades)
pred <- predict(glm_algas2, newdata, type = "response", se.fit = TRUE)
# Calcular intervalos de confianza
newdata$fit_link <- pred$fit
newdata$lwr_link <- pred$fit - 1.96 * pred$se.fit
newdata$upr_link <- pred$fit + 1.96 * pred$se.fit
# Convertir a escala de probabilidad (logística)
inv_logit <- function(x) exp(x) / (1 + exp(x))
newdata$fit <- inv_logit(newdata$fit_link)
newdata$lwr <- inv_logit(newdata$lwr_link)
newdata$upr <- inv_logit(newdata$upr_link)
# Promediar sobre el factor
pred_marginal <- aggregate(cbind(fit, lwr, upr) ~ Long, data = newdata, mean)
# Graficar relación promedio con IC
library(ggplot2)
ggplot(pred_marginal, aes(x = Long, y = fit)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "#0072B2", alpha = 0.2) +
geom_line(color = "#0072B2", linewidth = 1.3) +
labs(
x = "Talla individuo (cm)",
y = "Predicción probabilidad estructuras reproductoras",
) +
theme_minimal(base_size = 14)newdat <- expand.grid(
Long = seq(min(algas$Long), max(algas$Long), length = 100),
Sp = levels(algas$Sp)
)
# Predicciones de probabilidad
pred <- predict(glm_algas2, newdat, type = "response", se.fit = TRUE)
# Convertir de escala logit a probabilidad
newdat$fit <- plogis(pred$fit)
newdat$lwr <- plogis(pred$fit - 1.96 * pred$se.fit)
newdat$upr <- plogis(pred$fit + 1.96 * pred$se.fit)
ggplot(newdat, aes(x = Long, y = fit, color = Sp, fill = Sp)) +
geom_line(linewidth = 1.2) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.15, color = NA) +
labs(
x = "Talla individuo (cm)",
y = "Predicción probabilidad estructuras reproductoras",
) +
theme_minimal(base_size = 13)Comparaciones entre especies
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: glm(formula = Estado ~ Long + Sp, family = "binomial", data = algas)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Tom - Fra == 0 -0.5748 0.3453 -1.665 0.217810
Ver - Fra == 0 -1.8488 0.3693 -5.006 < 1e-04 ***
Ver - Tom == 0 -1.2740 0.3126 -4.076 0.000136 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)